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Abstract 

The origin of the spurious interface velocity in finite difference lattice Boltzmann 
models for liquid - vapor systems is related to the first order upwind scheme used 
to compute the space derivatives in the evolution equations. A correction force 
term is introduced to eliminate the spurious velocity. The correction term helps to 
recover sharp interfaces and sets the phase diagram close to the one derived using 
the Maxwell construction. medskip 

Keywords: Lattice Boltzmann; Liquid - Vapor Systems; Spurious Interface 
Velocity. 

1. General description of Finite Difference Lattice Boltzmann 
models 

Lattice Boltzmann (LB) models El El provide an alternative to current sim- 
ulation methods in computational fluid dynamics. These models are based on the 
physics at the mesoscopic scale, so that the macroscopic phenomena are recovered 
without solving the equations of continuous media mechanics. The starting point 
of any LB model is the Boltzmann equation |S] 

collisions 



dt m ) V dt 
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Here / = /(r, v, t) is the distribution function of fluid particles (supposed to be 
identical, with mass m), v is the particle velocity and F = F(r, t) is the local 
force acting on fluid particles. The collision term in the Boltzmann equation is 
usually linearized using the Bhatnagar-Gross-Krook (BGK) approximation j^j after 
introducing a relaxation time r: 



i>J ] = ~ [/(r,v,t) - /"*(r,v,t)] (2) 

collisions 



dt 



The equilibrium distribution function f eq = f eq (r, v, t) which appears in eq. (J2J 
is the Maxwell - Boltzmann distribution function: 

/-(r, v, t) = n(r, t) ( ^ ) ^ exp { - ^ • [v - u(r, ttf } (3) 

where ks is the Boltzmann constant, T is the absolute temperature of the system, 

n(v,t) - y /(r,v,t)dv (4) 

is the particle number density and 

u(r,t) = -^/v./(r,v,t)dv (5) 

is the local fluid velocity 0. We assume that the system is not too far from the 
equilibrium state and get [HUHj 

777 

V v /(r, v, t) ~ V v / e «(r, v, t) = - [ v - u(r, t) } / e «(r, v, t) (6) 

Recent investigations |X(JI 1111 IT2*] resulted in a general procedure to construct 
lattice Boltzmann models for single - component fluids. After discretization of the 
phase space |2 13 E] , the distribution functions are defined only in the nodes x 
of a discrete lattice C in the one - (ID), two - (2D) or three - dimensional (3-D) 
space, while the velocities are reduced to a discrete set {e^}, i = 0,1,., .TV. The 
elements of the ID, 2D and 3D velocity sets {e^} currently used in the LB literature 
[3J0IEI1EII are expressed using the propagation speed c = y'/s^T/xm, where x is a 
constant specific to each set. Following the discretization procedure, the Boltzmann 
equation (JJJ is replaced by the set of N equations 

d t fi(x, t) + e, • V/;(x, t) = - - lM*,t) - /f(x,t)] 

T 

' jf(x,i) [e t - u(x,t)] • F(x,t) 



Xc 2 " 



(i = 0, 1, ... AT) (7) 
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where the distribution functions /j(x, t) express the probability of finding at node 
x G C a particle having the velocity e^. The particle number density n and the local 
velocity u are now expressed as 



= n(x, t) = £ /<(*, i) = E /f ( x - *) 



(8) 



u = u(x, *) = ^ £ *) = ^)S ^"(x, *) (9) 

v ' y f=0 i=0 

while the equilibrium distribution functions in (J| are given as series expansion in 
the local velocity: 
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ej • u (ei • u) 2 u • u 

Xc 2 



2 X 2 c 4 



2 X c 2 



(10) 



In the one - dimensional D1Q3 model 4 , x = 1/3, M = 2 and the velocities ej 
are given by: 
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4/6 (« = 0) 
1/6 (i = 1, 2) 
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When using a square lattice in the 2D space (D2Q9 model Hp. v = 1/3, TV = 8 
and the velocities ej are given by 



r 



ej = < 



cos v 2 ; , sin v 2 ; 



(1 



(j-5)7T 
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while the weight factors are 



(i— 5)-7r 
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The Cartesian projections &i a = (e,) Q (i — 1,2,... , N; a — x, y, . . .) of the 
velocity vectors e< satisfy the relations 

} J w t e ia = 

i 

Wie iot ei0 = xc 2 (5 Q/ 3 

i 

i 

^2 w i e la e l! 3ei 1 ei(, = x 2 c 4 (S a pS lS + SpySs a + SajSps) (15) 

i 

The following sums are easily computed using the definition l|lUfl of the equilibrium 
distribution functions: 

£/r = n 

i 

i 

^ ei a eipf* q = n[xc 2 5 a p + u a U(3] 

i 

/] e-i<xeipe il fl q = nxc 2 [<5 Q/3 u 7 + 6^u a + <5 7Q m,3] (16) 

i 

Here w Q (a = x, y) are the Cartesian components of the local velocity u. 

The set of phase space discretized LB equations Q for the distribution func- 
tions fi = fi (x, t) may be solved numerically using an appropriate finite difference 
scheme defined on the lattice C. When using a scheme based on the characteristics 
line, the forward Euler difference is used to compute the time derivative while the 
first - order upwind scheme may be used for the space derivative ^3] as usually done 
in classical LB models 0] . These schemes give the following updating procedure 
for the distribution function ^2 ED 1X0] 

cSt 

fi(x,t + 6t) = /i(x,t) + — [/i(x,t) - /i(x- 8sei/c,t)} (17) 
— [/,(x,t) - /f (x,i)] + (x,i) [e, - u(x, i)] • F(x, t) 

T XC 

where St is the time step and Ss is the lattice spacing. 

2. Dimensionless equations and the force term 

Let In, tjt, tir, Tfl, cr and an the reference quantities for length, time, particle 
number density, temperature, speed and acceleration. In order to preserve the 
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form of the LB evolution equations, the reference quantities should satisfy the 
following relations 



tRCR 

Ir 

a R^R 
CR 



= 1 



1 (18) 



Since we will refer further to a van der Waals fluid, a natural choice for ur, Tr 
and cr may be the values corresponding to a mole of fluid at the critical point 
H3, riR = N A /V mc , T R = T c and c R = ^k B T c /m, respectively. Here N A is 
Avogadro's number, V mc is the molar volume at the critical point and T c is the 
critical temperature. With this choice of the reference quantities, the dimensionlcss 
propagation speed is 

c = (19) 

where T is the dimensionless temperature (note that T = 1 at the critical point). 
The other reference quantities may be derived from Eqs. (|18|l . if we choose the 
characteristic system size as reference length. 

The fact that the dimensionless propagation velocity c (|19l) should not neces- 
sarily equal 1 is a characteristics of the Finite Difference LB models where 
the propagation velocity is no longer related to the lattice spacing as in standard 
(classical) LB models Comparison of fluid velocity profiles at different tem- 

peratures becomes easier with the definition (|f 9J) associated to the expression Q of 
the local fluid velocity. 

The local force acting on a fluid particle is the sum of two terms: 

F(x, t) = F v + F a (20) 

The first term F v accounts for phase separation in the van der Waals fluid, while 
the second one F CT generates the surface tension at the interface between phases. 
The expression of F v and F CT is adopted from the literature [H] 

F^ = -V(- Pw + X c 2 p) (21) 
P 

F a = KV(V 2 p) (22) 

Here p w is the van der Waals pressure, x° 2 P is the dimensionless pressure of the ideal 
gas and k is a parameter that controls the surface tension. When using the reference 
pressure pr = mriRC 2 R , which is in accordance to Eqs. Q18[l. the dimensionless van 
der Waals equation of state reads 

Pw = ^--^P 2 (23) 
3 — p 8 

Note that the dimensionless form (|23() of the van der Waals equation of state is 
different from the form used by other authors 120] ■ 
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Figure 1: Plane interfaces in a 2D liquid - vapor system (the high density phase is black) 
Lattice size: 100 x 25 nodes; lattice spacing: 5s = 0.01; time spacing: 5t = 0.001 
temperature: T = 0.70; surface tension parameter: k = 0.00005. 
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Figure 2: Initial density profiles for two values of the dimensionless temperature. 



3. Stationary case: origin of the spurious velocity in the interface 
region 

Simulation results reported in this paper refer to the equilibrium state of flat inter- 
faces in liquid - vapor systems like the 2D one shown in Figure 1. Density, velocity 
and pressure profiles, as well as the phase diagram recovered in the stationary (equi- 
librium) state are found to be identical when using both D1Q3 and D2Q9 models. 
For this reason, we refer to the D1Q3 model on a ID lattice with N = 100 nodes 
(i.e., Ss — 0.01) when reporting the simulation results, as follows. 

For each value of the temperature, the initial density profile was set up to allow 
a transition region (Figure 2), while the values of the initial density of the liquid and 
vapor phases were computed using the Maxwell construction |17|. This procedure 
(which avoids large density gradients) helps the system to exhibit a stable behavior 
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Figure 3: Equilibrium density profiles recovered with the first - order upwind scheme 
(|17jl for two values of the dimensionless temperature and different values of the surface 
tension parameter k. 



until the equilibrium state is reached. The dimensionless relaxation time was set 
to t = 0.01 and we always performed 200,000 time steps (St — 0.001) using the 
upwind finite difference scheme Ijl7|) to ensure the equilibrium state of the fluid 
system. 

Simulations were done with three values of the parameter k which controls the 
surface tension. Figure 3 shows the resulting density profiles, for T = 0.90 and 
T = 0.80. As expected the interface width becomes larger as K increases. 

But there is still a transition region between the two phases for n = 0, when one 
would expect a sharp interface. Velocity profiles derived in the equilibrium state 
(Figure 4) show the existence of the spurious velocity (i.e., unphysical nonvanishing 
values of this quantity) in the interface region. The magnitude of the spurious 
velocity becomes larger when decreasing the temperature i.e., when the difference 
between the high and low density phases increases. However, the spurious velocity 
in the interface region reduces for larger values of the surface tension parameter k, 
when density gradients in the interface region become smaller. Also, the interface 
profiles of the van der Waals pressure (|23[) shown in Figure 5 become wider when the 
value of k is increased. For temperatures near the critical point and large values of 
k, the left and right interfaces overlap one another. This is especially seen in Figure 
5c, where the pressure plateaux of both the high and low density phase vanish for 
T = 0.90. 

The phase diagram of the liquid vapor system (Figure 6) is strongly affected by 
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Figure 3: (cont'd) Equilibrium density profiles recovered with the first - order upwind 
scheme (|17f) for two values of the dimensionless temperature and different values of the 
surface tension parameter k. 



S 



the value of the surface tension parameter k. The liquid - vapor system becomes 
more stable (i.e., lower temperature states may be reached) when using larger values 
of k which ensures smaller density gradients in the interface region. Moreover, when 
k is increased, the computed values of the densities (especially in the vapor phase) 
become closer to the theoretical values derived using the Maxwell construction. The 
phase diagram becomes flawed near the critical point for larger k, but this is an effect 
of the overlapping of the right and left interfaces. 

To account for the spurious interface velocity, we recall that the real LB equa- 
tions ^3] solved using the upwind finite difference scheme in the stationary case 
are, up to second order in the lattice spacing Ss 

e % ■ V/i(x, t) - ij;e l pe %1 dpd 1 f t {yi,t) = -— [f t (x,t) - f* q (x,t)} 

T 

{i = 0, 1, . . . AO (24) 

where ip — Ss / c. Consequently, the stationary mass and momentum equations 
recovered from Eqs. 124fl using the Chapman - Enskog procedure up to second 
order in the Knudsen number, are 

dp(pup) = ipdpd-f [xc 2 p8 l 3 1 + pupUj] (25) 
dpipuaup) + d a p w = xc 2 Tdf3 [pd a up + pdpu a ] (26) 

+ xc?ty [2d a dp{pu p ) + V 2 (pu a )] + P Kd a {V 2 p) 

From the mass equation (|25|l we get the expression of the spurious velocity which 
is present in the interface region where the fluid density is not constant: 

up = ^"^7 [x c2 /?Vy + pupu 7 ] ~ XC ^ — dpp (27) 

Thus, the spurious velocity is related to the numerical error introduced by the first 
order upwind finite difference scheme. Eq. (|27f) is in good agreement with simulation 
results (Figure 7) , especially when the system temperature is near the critical point 
or when the surface tension parameter k is large, i.e., when density gradients and 
the magnitude of the spurious velocity in the interface region are small. 
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Figure 4: Equilibrium velocity profiles recovered with the first - order upwind scheme 
l|17jl for two values of the dimensionless temperature and different values of the surface 
tension parameter k. 
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Figure 4: (cont'd) Equilibrium velocity profiles recovered with the first - order upwind 
scheme 1)17(1 for two values of the dimensionless temperature and different values of the 
surface tension parameter k. 
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Figure 5: Equilibrium pressure profiles recovered with the first - order upwind scheme 
()17() for two values of the dimensionless temperature and different values of the surface 
tension parameter k. 
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Figure 5: (cont'd) Equilibrium pressure profiles recovered with the first - order upwind 
scheme 1)17(1 for two values of the dimensionless temperature and different values of the 
surface tension parameter k. 
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Figure 6: Phase diagrams recovered with the first - order upwind scheme (|17|) for different 
values of the surface tension parameter k. 
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Figure 6: (cont'd) Phase diagrams recovered with the first - order upwind scheme Q17|) 
for different values of the surface tension parameter k. 




Figure 7: Velocity profile at T = 0.90: comparison between numerical results and 
Eq. (|27jl for different values of the surface tension parameter k. 
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Figure 7: (cont'd) Velocity profiles at T = 0.90: comparison between numerical results 
and Eq. (|27|) for different values of the surface tension parameter k. 
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4. Correction force term: reduction of the spurious interface ve- 
locity 

To eliminate the spurious term in the mass equation l|25p. we consider a supplemen- 
tary force term 



F?p = — - [-dp^eijAy) + (u^dp + upd^Aj] 



(28) 



where 



^7 = X<?&iP + ds(puyus) (29) 
After introduction of the correction force term (|28() , the LB equations (|24|l become 

1 



ej • V/i(x, t) - ipe^eijdpdjfi^t) 



[/i(x,t) - jf(x,*)] + 



• [e i/3 - up(x, t)][Fp + F» p ] (i = 0, 1, ... AT) 



(30) 



Introduction of the correction term cancels the spurious velocity term and allows 
the recovery of the correct mass equation in the stationary case. This is easy to see 
after calculating the sum 



J2 -^2 /f i e iP - U P) F h 

xs- 



— j — ff {&ip - up) [~dp{e il A 1 ) + (ujdp + upd 7 )A 7 } = (31) 
■ Xc P 

-ij)dpA 1 8p 1 = -ibxc 2 dp [dpp + d s (pupu s )] 

Moreover, the introduction of the correction force term does not alter the momentum 
equation since 

^ -^y — ft 9 e ia (eip - up) [-dp(e il A 1 ) + (u 7 dp + upd 1 )A 1 ] ~ 
i x ° p 

— ib [S a pu^ + Sp 7 u a + S ai Up] dpA 7 + ip5 a p [u-ydp + upd y ] A 7 + (32) 

ibSajdpA^ = 

In the ID case there is no need for Cartesian indices and the expression (|28() of 
the correction force term becomes 



F» = -^{e l -2u)[ X c^p + V 2 {pu 2 )} 
P 



(33) 



The following formula is used to compute the effect of the Laplace operator on a 
function / (e.g., p or pu 2 ) defined in the nodes x of the ID lattice: 

i=JV 

Y,Wif(x + eiSs/c,t) - /(x,f) CU) 



V 2 /(x,i) 



x(Ss) 



i=0 
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Figure 8: Profiles recovered using the correction force term (|28|) for T = 0.85 and 
K = 0. 



Simulations done with the correction force term l|28[l included in the LB evolution 
equations give sharp density profiles for k = (Figure 8a), as expected. The 
spurious velocity is no longer present in the interface region. A wiggly profile of the 
local velocity is still observed (Figure 8b), but the amplitude of the oscillations are 
several orders of magnitude smaller than the typical value of the spurious velocity 
arising in the interface region when using the uncorrected upwind scheme. Also 
the pressure profile recovered for k = is quite flat on the whole lattice, although 
wiggles of very small amplitude are still present (Figure 8c) . 

As previously observed during simulations with the uncorrected upwind scheme, 
increasing the surface tension always helps to stabilize the system for lower values 
of the temperature when the correction force term is considered in the LB evolution 
equations (Figure 9). However, the phase diagrams recovered using the corrected 
force term (Figure 9) are closer to theoretical results than the phase diagrams de- 
rived with the bare upwind scheme (Figure 6). For large values of the surface tension 
parameter k, the overlapping of the left and right interfaces is still present for tem- 
peratures close to the critical point, but the plateau is clearly defined in both high 
and low density phases for lower temperatures (Figure 10). The velocity profile 
recovered for T = 0.50 and k = 0.0001 is similar to the one shown in Figure 6b, 
but the amplitude of the oscillations is significantly reduced (2.0 x 10~ 14 ) because 
of the stabilization effect of the surface tension. 
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Figure 8: (cont'd) Profiles recovered using the correction force term (|28jl for T 
and k = 0. 
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Figure 9: Phase diagram recovered with the correction force term (|28|). 
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Figure 9: (cont'd) Phase diagram recovered with the correction force scheme 1)28(1 . 
5. Conclusions 

The first order upwind finite difference scheme used for the operator • V in the 
LB evolution equations Q generates the spurious velocity in the interface region 
of liquid - vapor systems. The magnitude of the spurious velocity is related to the 
density gradient through Eq. (|27|l . Consequently the spurious velocity is reduced 
by the surface tension. 

A correction force term (12811 is suggested to reduce the spurious velocity. Intro- 
duction of this term in the LB evolution equations allows to get sharp interfaces 
when the value of the surface tension parameter vanishes. Moreover, the phase 
diagram of the liquid - vapor system becomes closer to the theoretical one derived 
using the Maxwell construction. 
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Figure 10: (cont'd) Profiles recovered using the correction force term (|28|) for T 
and k = 0.0001. 
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